# Model one 
model{	
	for(i in 1:N){	# loop over observations 

	   # Utility Party 1
	   U11[i] <- beta[1,1] + delta[1] * (gamma[1] * UP[i,1] + (1-gamma[1]) * UC[i])  +   delta[6] *lr[i,1]  + delta[4] * pid[i,1] + X[i,] %*% psi[1,1,1:6] 
	   U21[i] <- beta[2,1] + delta[2] * (gamma[2] * UP[i,1] + (1-gamma[2]) * UC[i])  +   delta[7] *lr[i,1]  + delta[5] * pid[i,1] + delta[3] * a11[i]   + X[i,] %*% psi[2,1,1:6] 

	   # Utility Party 2
	   U12[i] <- beta[1,2] + delta[1] * (gamma[1] * UP[i,2] + (1-gamma[1]) * UC[i])  +   delta[6] *lr[i,2] + delta[4] * pid[i,2] + X[i,] %*% psi[1,2,1:6] 
	   U22[i] <- beta[2,2] + delta[2] * (gamma[2] * UP[i,2] + (1-gamma[2]) * UC[i])  +   delta[7] *lr[i,2] + delta[5] * pid[i,2] + delta[3] * a12[i] + X[i,] %*% psi[2,2,1:6] 

	   # Utility Party 3
	   U13[i] <- beta[1,3] + delta[1] * (gamma[1] * UP[i,3] + (1-gamma[1]) * UC[i])  +   delta[6] *lr[i,3]+   delta[6] *lr[i,1]+ delta[4] * pid[i,3] + X[i,] %*% psi[1,3,1:6] 
	   U23[i] <- beta[2,3] + delta[2] * (gamma[2] * UP[i,3] + (1-gamma[2]) * UC[i])  +   delta[7] *lr[i,3]+   delta[6] *lr[i,1]+ delta[5] * pid[i,3] + delta[3] * a13[i] + X[i,] %*% psi[2,3,1:6] 


	   # Utility 4
	   U14[i] <- 0
	   U24[i] <- 0

	    # Normalize

	    P1[i,1] <- exp(U11[i]) / (exp(U11[i]) + exp(U12[i]) + exp(U13[i]) + exp(U14[i]))  
	    P1[i,2] <- exp(U12[i]) / (exp(U11[i]) + exp(U12[i]) + exp(U13[i]) + exp(U14[i]))  
	    P1[i,3] <- exp(U13[i]) / (exp(U11[i]) + exp(U12[i]) + exp(U13[i]) + exp(U14[i]))  
	    P1[i,4] <- exp(U14[i]) / (exp(U11[i]) + exp(U12[i]) + exp(U13[i]) + exp(U14[i]))  

	    P2[i,1] <- exp(U21[i]) / (exp(U21[i]) + exp(U22[i]) + exp(U23[i]) + exp(U24[i]))  
	    P2[i,2] <- exp(U22[i]) / (exp(U21[i]) + exp(U22[i]) + exp(U23[i]) + exp(U24[i]))  
	    P2[i,3] <- exp(U23[i]) / (exp(U21[i]) + exp(U22[i]) + exp(U23[i]) + exp(U24[i]))  
	    P2[i,4] <- exp(U24[i]) / (exp(U21[i]) + exp(U22[i]) + exp(U23[i]) + exp(U24[i]))  

	   # Probailities
	   P[i,1]   <-  P1[i,1] * P2[i,1]
	   P[i,2]   <-  P1[i,1] * P2[i,2]
	   P[i,3]   <-  P1[i,1] * P2[i,3]
	   P[i,4]   <-  P1[i,1] * P2[i,4]

	   P[i,5]   <-  P1[i,2] * P2[i,1]
	   P[i,6]   <-  P1[i,2] * P2[i,2]
	   P[i,7]   <-  P1[i,2] * P2[i,3]
	   P[i,8]   <-  P1[i,2] * P2[i,4]

	   P[i,9]   <-  P1[i,3] * P2[i,1]
	   P[i,10] <-  P1[i,3] * P2[i,2]
	   P[i,11] <-  P1[i,3] * P2[i,3]
	   P[i,12] <-  P1[i,3] * P2[i,4]

	   P[i,13] <-  P1[i,4] * P2[i,1]
	   P[i,14] <-  P1[i,4] * P2[i,2]
	   P[i,15] <-  P1[i,4] * P2[i,3]
	   P[i,16] <-  P1[i,4] * P2[i,4]


	   # Distribution
	   Y[i] ~ dcat(P[i,1:16])

       }

       	 	  
       gamma[1] ~ dunif(0,1)
       gamma[2] ~ dunif(0,1)
	


	# Other Prior
	delta[1] ~ dnorm(0,.001)
	delta[2] ~ dnorm(0,.001)	
	delta[3] ~ dnorm(0,.001)
	delta[4] ~ dnorm(0,.001)	
	delta[5] ~ dnorm(0,.001)
	delta[6] ~ dnorm(0,.001)	
	delta[7] ~ dnorm(0,.001)


	beta[1,1] ~ dnorm(0,0.001)
	beta[1,2] ~ dnorm(0,0.001)
	beta[1,3] ~ dnorm(0,0.001)
	beta[2,1] ~ dnorm(0,0.001)
	beta[2,2] ~ dnorm(0,0.001)
	beta[2,3] ~ dnorm(0,0.001)

	# psi
	for(k in 1:6){
	psi[1,1,k]  ~ dnorm(0,.001)	
	psi[1,2,k]  ~ dnorm(0,.001)	
	psi[1,3,k]  ~ dnorm(0,.001)
	psi[2,1,k]  ~ dnorm(0,.001)	
	psi[2,2,k]  ~ dnorm(0,.001)	
	psi[2,3,k]  ~ dnorm(0,.001)		
	}

	# Sample Difference
	diff <- gamma[1] - gamma[2]

}
